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The thesis presents development of the nonlinear 


differential equations both in stress function-vertical 
displacement (0-w) and displacements in X— , Y — and In- 
direction (u-v-w) forms for a hyperbolic paraboloid. As it 
is found difficult to solve them by analytical means, 
recourse has been taken to finite element method. , The 
shell has been discretized into thin shallow doubly curved 
rectangular elements. 


The stiffness matrix for a thin shallow doubly 
curved rectangular element with five degrees of freedom 
at each node has been generated by using numerical 
integration technique with 6-point Gauss-Legendre quadrature. 
This is independently checked by using a technique suggested 
by Kabaila. 



The linear bending analysis using the above stiffness 
matrix has been' applied to hyperbolic paraboloid shells 

I" 

bounded by A4;s characteristics under uniform vertical load 
for the following boundary conditions: (i) Simply supported 
and (ii) Clamped. The deflections, in-plane shears, bending 
moments etc., have been plotted for the various boundary 
conditions. The computer program which forms the most 
important tool in the finite element analysis has been 
developed independently. The convergence of results has 
been illustrated for 4x4, 6x6 and 8x8 mesh sizes of a 
hyperbolic paraboloidal shell with fixed edges carrying a 
uniform vertical load. 

The nonlinear bending analysis for the case of 
large deflection (i.e. geometrically nonlinear) has been 
formulated and applied to a hyperbolic paraboloid shell 
with fixed edges under uniform vertical load using a 
4x4 mesh size. The study has been confined to this 
particular mesh size in view of time and memory constraints 
on the computer. Newton-Raphson iteration has been adopted 
to trace the loa.d-deflection curve at the centre of the 
shell. The tangent stiffness matrix and the geometric 
stiffness matrix have been generated by numerical integra- 
tion technique by writing a special computer program. The 
buckling load for the shell has been indicated from the 
load-deflection behaviour and compared with the existing 
experimental results.. 


Four appendices have been provided to cover 
(i) numerical integration technique using Gauss-legendre 
quadrature, (ii) computer integration of congruently 
transformed matrices due to Kabaila, (iii) stiffness 
matrix of straight beam with respect to eccentric axes 
and (iv) skewed boundary conditions. 


CHAPTER 1 
INTRODUCTION 

1.1 GENERAL 

Shallow doubly curved translational shells are 
efficient and aesthetically attractive structures and as 
such are being used extensively in roof construction. Due 
to the limited applicability of the membrane or 'momentless' 
theory, the bending theory of shells has received a great 
deal of attention in recent years (1,2, 5, 4,5.6) *. 

The* shallow shell theory' of Margeurre (1) and 
Vlasov (2) is often used as it is reasonably accurate for 
the shell dimensions used in practice and since it is 
considerably simplified in comparison with more acourate 
shell theories. The fundamental approximation of the shallo 
shell theory is that quadratic .terms involving the slopes oi 
the shell middle surface are negligible compared to unity. 
While there is certain uncertainty regarding the magnitudes 
of the errors introduced by the approximation, a common 
rule-of-thumb is that if the maximum rise/span ratio of 
the shell is less than one-fifth, the theory produces 
results sufficiently for practical purposes (83). 


* Numbers in parentheses refer to references listed in 
the list of references. 
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The shallow shell theory is normally formulated in 
terms of either a stress function and normal displacement 
(0-w) or in terms of middle- surface displacements u,v and 
w. The differential equations for both the above formula- 
tions have been derived for the large deflection case -along 
with the boundary conditions. In the case of hyperbolic 
paraboloid, the harmonic functions do not suit because of 
the asymmetric relation between w and 0 (7). Thus, despite 
the simplifying assumptions which are made, problem is 
still difficult one mathematically with the result that 4 


comparatively few analytical solutions have been obtained * 

A recent study (8) has quantitatively demonstrated- the 
sensitivity of shells to boundary conditions thus emphasi- , 
zing the importance of developing methods of analysis which; 


are not restricted to special boundary conditions, parti- 
cularly those specified by the simple support case. 
Evidently the most promising approach to the problem is a 
numericaljone. Among the numerical procedures a distinc- 
tion is often made between procedures which Eire regarded 


as mathematical approximations, such as finite differences 
and variational methods, and those which are physical 
approximations, such as various discrete element systems- 


Das Gupta (9), Soare (10), Mirza (11), Russell and 
Gerstle (12,13) and others have presented finite difference 
solutions for the linear analysis of shallow shells 
particularly as related to various hyperbolic paraboloids? .. 


as bounded by characteristics. Abu-Sitta has applied 
finite difference methods to elliptic paraboloids (14). 

In the application of the variational methods to 
the analysis of shallow shells (7), the problem of 
selection of approximating functions restricts the 
versatility of the method since a definite choice of such 
functions applies only to specific boundary conditions. 

In addition, the treatment of non-classical boundary 
conditions presents serious difficulties. 

Among discrete element procedures which have been 
applied to shell problems are methods using the lumped 
parameter model (15,16), the framework or lattice model 
(17) and the finite element technique. While it is 
justified to classify the first two methods as physical 
idealizations*, in some applications ** the finite element 
approaches to plane stress (18) and plate bending (19,20) 
have met with notable success. 


* The lumped parameter model may be regarded as a 

physical interpretation of finite difference approxi- 
mations since the governing equations of the model 
are equivalent to the consistent difference equations 
of the problem. However, the model facilitates the 
formulation of boundary" conditions without using 
fictitious grid points. 

** The applications referred to are those in which the 
shape of the structure is not idealized. When*, for 
example, flat elements are used to approximate a 
curved surface, a physical idealizatiOh is also 
involved . 


4 


A natural extension of the scope of the finite 
element method is to shell problems. Most of the early 
attempts to extend the method employed assemblages of 
flat elements to approximate the curved surface of the 
shell (21). Such an idealization may not be entirely 
satisfactory, however, since errors are introduced which 
are distinct from those involved in assuming the form of 
displacements in the structure (or element). Thus no 
direct relationship with the Rayleigh-Ritz procedure is 
apparent. Jones and Strome (22) demonstrated the need 
for the use of curved shell elements when solutions were 
sought for shells for revolution. Curved elements have 
been used for axisymmetric shells (23) and, for the cases 
that have been reported, more accurate results have been 
obtained than with flat elements. In the case of trans- 
lational shells, the use of the doubly curved element 
gives more accurate result than that of flat element even 
for coarse mesh size. Since hoop stiffness present in 
shell of revolution problems is frequently not present in 
many translational shells, Schnobrich et al (21) are- of the 
view that flat elements may be used for their analysis. 

This means that a curved surface is idealized as a folded 
plate. To derive a stiffness matrix for a completely 
general doubly curved shell element would be rather 
difficult because of the complicated differential geometry 
that would be involved and the operations on large size 
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matrices that would he required in the derivation. Because 
of this difficulty, Schnohrich et al (21) have gone in for 
flat triangular elements, in which there is no coupling 
between membrane stiffness and plate bending stiffness. 

This results in the simplicity of the derivation of the 
stiffness matrix making it possible to use higher order 
displacement functions to represent in- plane and bending 
displacements. If the number of elements which are Used 
to represent the shell remains unchanged, using high order 
displacement functions for each element can efficiently 
improve the accuracy of the solution and the continuity of 
stresses between the elements. However, the corresponding 
increase in the number of generalized displacements in each 
element causes considerable increase in the number of 
simultaneous equations which have to be solved. This 
results in an increase in the band width and the computer 
storage capacity required. The computational effort 
becomes considerable and puts a serious limit on the 
number of elements which can be used. In this case, the 
advantages that can be gained by using higher order dis- 
placement functions tend to be destroyed by the poor geometri 
approximation growing out of the restriction on the number of 
elements that can be used. Thus, when effort is expanded on 
developing or choosing a flat shell element, it is important 
to consider the best balance between the displacement 
approximation and the geometric approximation. 
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Connor and Brebbia (24), Dhatt (25) and Pecknold 
and Schnobrich (26) have developed doubly curved finite 
elements using shallow shell theory and comparatively 

; 

low order displacement functions. Use of curved elements 
for translational shells of constant curvature involves no 
geometric approximation. 

1.2 REVIEW OP LITERATURE 

As Kovach (27) puts it, it is necessary to analyse 
for the nonlinear phenomena that abound so widely in 
nature for clear understanding of the surrounding environ- 
ment. The complexities of nonlinear analyses have created 
a 'non-linear barrier* which, only in very recent times, [ 
has begun to be penetrated through the efficient use of j 

modern high speed computers. In structural mechanics, the 

'■ . . ■ - ' ■ ' ' ■ . . I 

unparalled success of the now well-known finite element 

| 

method gives hope that the nonlinear barrier can, at last, 
be broken. Depending on the sources of nonlinearities, [ 

nonlinear problems can be divided into three categories. | 
In brief, these categories ahe problems involving material 
nonlinearity alone, problems involving geometric nonlineari- 
ty alone, and problems involving both geometric and material 
nonlinearities. 

1.2.1 Nonlinearities 

Material nonlinearity is the easiest to visualise. 

It encompasses problems in which stresses are not linearly 
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proportional to the strains, hut in which only small dis- 
placements and small strains are considered. Displacements 
refer to the changes in the overall geometry of the body 
whereas strains are related to internal deformations. The 
word 'small* usually implies infinitesimal changes in the 
geometry of the body. An example is elastic-plastic 
analysis of structures. Linear strain-displacement relations 
are used. 

Geometric nonlinearity includes problems arising 

■■ ■ ■ '■ ' '■ ' ■ ■ . .. '' ' . ■ [ 

both from nonlinear strain-displacement relations and 
finite changes in geometry. Here linear stress-strain 
relations are used. In other words, this category encom- 
passes large strains and large displacements. An important 
subclass of geometrically nonlinear problems is the case of 

small or infinitesimal strains and large or finite displace- ; 

j 

ments. An example is the elastic post-buckling behaviour 
of structures. j 

Finally, the third and most general category of j 

nonlinear problems is the combination of the first two I 

categories. It involves nonlinear constitutive behaviour ! 
as well as large strains and finite displacements. The 
deformation of a rubber like material is an example of the 
third category. 

In large span shells. It becomes necessary to cut 
down the weight of the structure by making the shell 
thin. This necessitates the analysis of shells for large 
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deflection taking the strains as small. Similarly the air- 
craft and ship structures are designed. The above come under , 
the problems of geometric nonlinearity which is now reviewed. 

1.2.2 Geometrically Nonlinear Analysis 

Various surveys of the present topic have appeared 
in recent years. Martin (28), gave a thorough review of all 
work in geometrically nonlinear finite element analyes 
upthrough 1969, the present review begins from there onwards. 
Oden (29) gave a correspondingly thorough assessment 
that covered a wide range of nonlinearities, including 
material nonlinearities and many special formsof higher- 
order nonlinearities (e.g., large strains.) which are excluded 
from the present review. Stricklin and his associates at 
the Sandia Corporation have contributed a great deal to the 
present topic in recent years ( 30 ). 

The formulation of element force-displacement 
equations for geometrically nonlinear analysis can be 
approached from many standpoints. The attention is restric- 
ted to the concept of stationary potential energy and to a 
Lagrangian discription of behaviour, i.e., displaced points 
are referred to a fixed set of axes. 

The typical governing equation of geometric nonlinear 
problem in indicial form is 

P i = N ij + ^ijk A ;j A k + 3 N ijkj|_ A j A k A * (1.1a) ^ 

or, in matrix form 
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m - [K] {A} + [N 1 (a)]{A} + [N 2 (A)]{A} (1.1b) 

■where [Z] is the usual linear stiffness matrix, 
and [N.^] an ^ are the first- and second- order 

'geometric' or 'incremental' stiffness matrices 
respectively. 

The elements of [N-^] are linear functions of the 
displacements while the elements of [Ng] are quadratic 
functions of the displacements. Thus, the above formula- 
tion ends up in the solution of nonlinear equations. The 
parameters to be solved for differ from one physical 
circumstance to another and they do not pertain exclusively 
to a tracing of the load-displacement response along a 
continuous curve, starting from the linear regime at small 
loads onto significant nonlinearities at higher loads. 
According to Gallagher and Mau (56), three areas of 
interest with respect to algorithmic tools are perceived: 

(l) general nonlinear analysis with reference to bifurcation 

f"£) 

phenomena; 'calculation of bifurcation or 'branching' points 

3 

and determination of the load -displacement response 

along a post-buckling path. 

1.2.3 General Nonlinear Analysis 

General nonlinear analysis using procedures which 
operate directly on the potential energy have two advantages 
Firstly, the representation of the system can be conveniently 
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established without recourse to large-scale matrix operations, 
since energy is a scalar. Secondly, it had been hoped that 
the calculation of a stationary value for a more general 
nonlinear function might not be significantly more difficult 
than that for a quadratic function. Bogner, Pox and Schmit 
(33) used function minimization procedures in nonlinear 
finite element analysis. It has been claimed that the 
energy search method may be used effectively for the geometric 
nonlinear analysis. Unfortunately, other researchers (for 
example, Oden and Key (34)) have reported disappointing 
results with the energy search method. Computational 
experience has not proved conclusively that the advantages 
of directly operating on the potential energy are indeed 
manifested in application. A view of the function minimiza- 
tion schemes is given by Broyden (35), but without reference 
to physical applications. 

Methods which attack directly the solution of the 
governing equilibrium equations are overwhelmingly the most 
popular in application to geometrically nonlinear analysis. 
Many different specific forms are available, including 
direct iteration, Newton-Raphson, incremental and initial 
value methods. In discussing these methods, it must be 
borne in mind that in tracing a load-displacement response 
through a significant range of nonlinear behaviour, it is 
essential to deal with intervals of loading, the solution 
for an interval being employed as the starting data for the 
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second interval. Numerical calculation shows that it is 
not possible, in general, to proceed from the unloaded to 
final state in a single step due to numerical instabilities 
in the solution. 

1.2.4 Solution by Direct Iteration 

Direct iteration is seemingly the most attractive 
approach to geometrically nonlinear analysis since it 
requires only the inversion (or solution) of the linear 
portion of the governing equations. The iterative 
sequence for the solution continues until the difference 
between two successive solutions converges to a specified 
tolerance. Convergence difficulties are encountered when 
the nonlinearities are severe. Such difficulties are 
often manifested by continued iteration in a loop about 
the convergent solution. In such cases an improved 
procedure is to employ a higher order iterative scheme 
as described by (36). The approach, in any case, suffers 
from inefficiency in comparison with certain alternatives. 
To improve this situation, Roberts and Ashwell (37) 
form the right hand side with use of an estimated mid- 
increment displacement field. 

1.2.5 Newton-Raphson Method 

Direct iteration can be interpreted as an approxi- 
mate form of the Newt on-Raphson method. To explain this 
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approach, all the nonlinear terms of the potential energy- 
are indicated hy Stricklin, et al (30). have shown 

that this method is of the form 

(tK] + rr”f »- )6 A = -W A °> -{ life + { p) (1.2) 

i 3 A j 9A 

where is formed on the "basis of {A 0 > and (6 A} 
symbolises the displacement interval covered by the 
iterative process. The solution for {6 A} is used to 
improve { A^ } and the process is continued until { 6 A } 
is sufficiently close to zero. From the efficiency stand- 
point, a major difficulty in the Rewton-Raphson approach 
is the need to reform the matrix on the left hand side and 
to invert it or solve the equations in each iteration. 

Also the step size for each solution is very critical 
since the method will not converge unless it is sufficien- 
tly close to the solution point. The method possesses 
second-order convergence properties. The Newt on-Raphson 
scheme may be modified to eliminate or minimize the 
repeated matrix inversion. 

1.2.6 Incremental Method 

The incremental method is probably the most 
widely-used of the procedures addressed directly to the 
solution of the equilibrium equations. It has been given 
a theoretical footing in the variational context by many 
authors, e.g. Hofmeister, et al (38 )» Fian and Tong (39) 
and Fellippa (40). In this, a tangent stiffness is 
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constructed for the interval in question and the solution 
for the interval begins with the solution of the prior 
interval as a base. Stricklin, et al (30) and Zienkiewicz 
and Nayak (41) also discussed about the same. Dupuis, 
Pfaffinger and Marcal (42) have made effective use of 
incremental method in shell applications, where initial 
imperfections are of considerable importance. 

1.2.7 Perturbation Method 

Walker (43) and Morin (44) have applied the 
perturbation approach to nonlinear analysis. Their 
results demonstrate the desirability of a corrector step, 
e.g., N ewton-Raphs on iteration. Although perturbation 
concepts show no particular advantage in this type of 
analysis, they are especially important in post-buckling 
analysis. 

1.2.8 Buckling Analysis 

The bifurcation (buckling) analysis is founded 
on the assumption that the distribution of internal 
forces in the structure, due to applied loads which cause 
bifurcation, can be computed using linear analysis. This 
procedure is, therefore, termed linear stability analysis. 
The second variation of the potential energy is zero for 
bifurcation (neutral equilibrium). Thus conservative 
loading is implied. Only a limited amount of work has been 
reported on problems of finite element elastic instability 
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analysis for non- conservative loadings (45). Contributions 
to finite element eigenvalue analysis on the iterative side 
have recently been made by Bronlund (46), Whetstone and 
Jones (47), Rosen and Rubenstein (48) and Dong et al (49). 
Method of conjugate gradients (50) has been applied to 
eigenvalue determination. 

The perturbation method and its variants constitute, 
almost exclusively, the approach to tracing load-displacement 
behaviour along the post-bifurcation branch. 

Geometrically nonlinear finite element analysis is 
on a sound theoretical footing. The basic relationships 
are generally accepted, many different plate and shell 
elements are capable of formulation for such analysis and 
solution algorithms abound. Computational efficiency and 
reliability are questionable areas, however. Dispropor- 
tionately small efforts have been devoted to reducing 
solution costs. The work of Stricklin and his colleagues 
( 30 ) is among the exceptions. 

The above review limits attention to procedures in 
which updating of structural geometry is not a critical 
aspect of the problems solved by these procedures. The 
work of Murray and Wilson (53) and allied developments are 
based upon the updating of geometry. Such updating is 
of course useful in enhancing the computational efficiency 
for moderate nonlinearities and is essential when the 
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nonlinearity is severe. Reference (30) studies the 
limitations of the subject procedures in this regard . 

Perhaps the most important class of practical 
problems for which geometric nonlinearity is important is 
that concerned with elastic instability. A hierarchy of 
finite element approaches to this problem is taking form. 

At one end of the scale is the simple approach of linear 
stability analysis. This approach gives an overestimate 
to the strength of imperfection sensitive structures. 
Almroth and Brogan (54) have very recently shown that 
equally serious underestimates are given for cases where 
sensitivity to initial imperfection is not a factor in the 
determination of the limit load. At the other end of 
the scale are problems governed by higher- order terms in 
the geometric nonlinearities and material nonlinearities. 
The papers by Oden (55) and Stricklin, Haisler and Von 
Riesemann (56) define well the current status of finite 
element approaches to their solution. 

1.2.9 Use of Hybrid Elements 

In 1964* Pian (72) suggested a method which was 
later discussed more rigorously by Pian and Tong (73) in 
which the inter element boundary displacement compatibility 
can be satisfied rather easily. In this approach, one 
assumes an equilibrium stress field in the interior of the 
element, and a displacement field at the boundary of the 
element, which inherently satisfies the interelement 
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compatibility condition. Pian's approach has since been 
employed in several small-deflection elastic problems such 
as plate bending, shell problems and analysis of multi- 
layered plates and shells where transverse shear effect 
is important. They clearly demonstrated the versatility 
of the hybrid stress model for deriving stiffness matrices 
of elements of arbitrary geometry. 

Using the above approach, some satisfactory results 
were obtained for incremental analysis of large deflection 
problems by Pirotin (74). As discussed by Pian (75), these 
methods cannot be considered as being based on the modified 
complementary energy principle which is the basis of the 
hybrid stress model and hence cannot be considered as 
consistent hybrid stress finite element methods. 

\ 

Since, as discussed by Pian (75), the hybrid assumed 
stress finite element model has proved to be more versatile 
structural tool, it is necessary to present a consistent 
variational formulation of the hybrid stress model for the | 
incremental analysis of large deflection problems. The 
concept of initial stress is to be employed, wherein, j 

during any given step in which the displacements, stresses j 

and external loads undergo increments, the state at the 
beginning of the step is considered as one of initial stress. 

A 'eheck 1 on the equilibrium in the state prior to the 
addition of further load increment is to be included. i 
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Since the hybrid stress model (76) has proven to be 
a valuable tool in the linear analysis of complex problems 
such as sandwich plates and shells and problems with 
singularities, it can be expected to be of equal value in 
large deflection problems. 

1-3 BASIC EQUATIONS 

The objective is to develop the nonlinear differen- 
tial equations for thin shallow shells having large deflection. 
The equations are specialized for a hyperbolic paraboloid 
bounded by elastic edge beams along the characteristics. 

1.3*3 General 

Medium bending of shell is defined (57) as that when 
the maximum deflection is of the same order as thickness or 
larger, but is small in comparison with other linear 
dimensions of the shell. As the shell is shallow, the 
effects of tangential displacements on rotations and 
change of curvatures is neglected. The assumptions of the 
shallow shell theory are: 

(A) The shell or plate is thin so that shear deformation 
is neglected. 

(B) Kirchhoff's hypothesis is accepted, i.e., normal 
stresses perpendicular to the middle surface are neglected 
and normal to the middle surface does not undergo any de- 
formation. 

(C) Medium bending of shpll is considered. 
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(D) The squares of elongations and shears and the 
products of elongations and rotations are neglected. The 
rotation about Z-axis is neglected. The squares of 
rotations about and which are of the order of 
elongations and shear are retained. 

(E) A shell is considered shallow when the rise to 
span (smaller) ratio is less than 1/5. Also L 2 /IL^R^ <<1 
where L is the bigger span. 

(E) In the case of shallow shells, Lome's first approxi- 
mation is used, i.e., (l+k^Z) * (l+k^z ) « 1. 

The quadratic terms involving the slopes of the shell 
middle surface are negligible compared to unity. If the 
shell middle surface is a second order surface, this 
assumption leads directly to the approximation that the 
curvatures of the middle surface are constant. An additional 
consequence is that the geometry of the surface is, in effect, 
approximated by that of its projection on the horizontal 
plane. Thus, if the parametric representation of the shell 
middle surface is given in the form 
X = X 

Y = Y (1.5) 

Z = Z(X,Y ) 

where (X,Y,Z) are cartesian coordinates, the approximation 
is made that the (X,Y) coordinate lines are orthogonal on 
the middle surface. It is in this sense that the term 
'orthogonal coordinates' will be used throughout the 
remainder of this study. 
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(G) Gauss- Codazzi equations hold good for the continuity 
of the surface. 


1.3.2 Strain-Displacement Relations in Orthogonal 
Curvilinear Coordinates 


let a. and k. . he the quantities which characterize 
the middle surface of the shell before deformation, referred 
to orthogonal curvilinear coordinates ^ During 
manufacturing and erecting of the shell, initial deforma- 
tions occur invariably. The irregularities may appear 
prior to the application of the load or owing to non- 
uniform temperature distribution in the shell. Initial 
stresses may also appear. These lead to deviation from 
orginal surface. It is assumed that the new surface is 
also shallow, so that the effect of tangential displace- 
ments on the initial rotations and the changes of curvature 
may be neglected. The superscript 'O' refers to initial 
deformations. 
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The elongations and rotations due to surface load may 
he determined by the following formulae (57) 
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The projections of the total displacements, the 
corresponding elongations and the changes of curvature will 
he summation of the quantities owing to initial deformations 
and the corresponding quantities due to load. 


1.3»3 Derivation of Equilibrium Equations in Orthogonal 
Curvilinear Coordinates 

If the stresses appearing in the shell during the 
manufacture and assembly are removed, but the initial 
deformation remains, then the irregularities are residual 
and that they appear as single normal displacement using 


The equations for the state of equilibrium of the 
shell under the normal load are derived using Reissner's 
variational theorem . The elongations at any point ’Z* from 
the middle surface are given by the following: 
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The expressions for the components of stresses at any point 
' z ' from the middle surface of the shell are as below: 
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where are s ^ ress resultants and M^,M ,M are 

stress couples. 

The variational equation may be written as 
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Substituting for the quantities in the above, the variation 
is carried out with respect to the unknowns u^, u^, w, 

N 2* N 12’ M l’ M 2 and M 12* Since 6u l’ du 2 * 6w ' 6N l’ 61 V m i2’ 
6M-^» and SM.^ are arbitrary, the coefficients of each 

must be zero. The above gives three equilibrium equations 

in u ^»u^ and w directions and six force-displacement 

relations along with proper boundary conditions. 
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The force-displacement relations are the standard ones 
and hence are not repeated here. 

I. 3.4 Reduction to a System of Two Nonlinear Equations 

A force function 0 is introduced which will satisfy 
the first two equilibrium equations of (1.9). 
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After introducing equations (1.10) in the first 
two equations of (1.9) » the left hand side of the two 
equations of (1.9) becomes 
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which can he neglected for shallow shell with non-zero 
Gaussian curvature and is zero for shells of zero Gaussian 
curvature. 

Also, 


EL + N = V 2 0 
12 


( 1 . 12 ) 


where operator v is defined as 
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A system of two nonlinear equations can he obtained 
in 0 and w from Gauss equation of continuity of def orma.tiori 
and third equilibrium equation (1.9). Gauss equation of 
continuity of deformation becomes 
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where k?. = k. . +x — 

ia 


The third equilibrium equation is reduced to an 


equation in w and 0 as follows: 


D ^ 2 V 2 w + IT-^ ^ + » 2 I + 2» 12 X 12 - p, 


( 1 . 16 ) 


where K 1} » k.. + X ± . + X y 


In the further analysis, initial deformations are 
not considered. In that case, the more general one, when the 
coordinate lines are not the lines of principal curvature, 
takes intoaccount the twisting of the surface and the two 
fundamental equations become 




28 


3 2 W \2 


Eh V * ^ “ { ( 3~5 ± 3^' 


9 2 w 

3 5?" 


3 2 w 

Ti 2 * 

2 


- k 


3 2 


w 


0 


D V 4 w - (k ? ? 


3 2 


£ 


a 2 d> 


+ v _ ot 

2 " 5 X l 3 5| 2t 3 


( 4 # 


3_f 0_ 


9 2 


w + 9 2 3 2 w 


"8 5 2 ^ 3S -L 3^2 3 Ej35 2 ' 85" H T 


i sir - k 2 


9 2 TJ- CS 2. 


y + 2t -} 

9 Cl3^ 


(1.17) 


) - P =0 
; N 


If the equation of the middle surface of a shallow shell has 
been assigned as Z = Z(5^, then it is possible to 

calculate approximate curvatures by using the approximate 
formulas 
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When the coordinate lines become the lines of principal 
curvature, the twist term disappears and the fundamental 
differential equations of nonlinear shell theory become 
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1.3.5 Specialization to the Doubly-Curved Shallow Shells 
Referred to Cartesian Coordinate System 

As the shell is shallow, the radii of curvature are 
constant and Lame's parameters are equal to Unity. The 
middle surface of a d oubly- curved shell is defined by 
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1.3*6 Stra.in— Displacement Relations 
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1.3*7 Equilibrium Equations 
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Force-displacement relations are of standard form and 
hence are not repeated here. 


1.3*8 Reduction to Two Nonlinear Equations in 0-w 

A force function 0 is introduced so that it satisfies 
the first two equilibrium equations of (1.22) 
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The equations can he reduced to a set of two nonlinear 


equations in 0 an d w. 
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1.5*9 Single Nonlinear Differential Equation in w for the 
Case of General loading 

When the coordinate lines are along the principal 
curvatures, the differential equation becomes 
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1.5.10 u-v-w Formulation of the Differential Equations 
When the coordinates are along the principal 


curvatures, 
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1.3*11 Elastic Edge Beam Boundary Conditions 


Referring to Figure 1.4, the following boundary 
conditions are deduced for edge beam parallel to Y-axis 
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1*3 »12 Nonlinear Differential Equations Specialized with. 

respect to Hyperbolic Paraboloid Shell Bounded by 
Characteristics 


From Figure 1.5, the equation for the shell middle 

f 

surface is given Z = ^ XT. Equations (1.24) become (58) 
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where k, = t 
1 ab 

1.4 OBJECT AND SCOPE 

The objective of this study is to apply the finite 
element method to the bending analysis of thin shallow 
hyperbolic paraboloid shells bounded by characteristics to 
a wide variety of boundary conditions including skew ones. 
Doubly-curved rectangular elements of negative Gaussian 
curvature with five degrees of freedom at each node have 
been used in the analysis. The analysis has been extended 
to the study of geometrically nonlinear behaviour of the 
above shell with fixed edges. The same may be applied to 
hyperbolic paraboloid shells with elastic edge beams. 
Attention has been confined to static analysis only 
assuming a homogeneous and isotropic shell material. 
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Normally, torsional rigidity (9) and the inplane bending 
stiffness of edge beams (12) are neglected as well as their 
eccentricity, if any, with respect to the shell middle 
surface. Inclusion of these fa,ctors presents no difficulty 
in the present analysis. 

The whole stiffness matrix is generated by numerical 
integration method using six point Gauss-Legendre quadrature 
since the stiffness matrix for the nonlinear case also is 
to be evaluated numerically. Attempt is also made to use 
the stiffness for shallow doubly-curved shell as developed 
by Connor and Brebbia (24) and also by Schnobrich and Pecknold 
(26) to verify the results. Since there weresome errors in 
the stiffnesses as seen from the computer results, the 
attempt is dropped finally as any corrective measure involves 
the derivation of the entire stiffness matrix which takes 
lot of time. 

The important step in the finite element method is 
the selection of appropriate displacement function. The 
two important displacement functions being used in plate 
bending analysis are discussed and the twelve-parameter 
polynomial expression for w used by Melosh (20 ) , 

Zienkifiwicz (19) and many others in plate bending analysis 
has been adapted here. 
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A variety of numerical examples is presented. 

In order to establish the validity of the method, 
comparisons are made with analytical solutions for 
hyperbolic paraboloids bounded by characteristics, 
both simply-supported and fixed, along the edges. Limited 
comparisons are made with other numerical results to 
further substantiate the reliability of the analysis . 

The dimensions of the earlier examples are retained for 
comparison of results. 

The programming part of the problem assumes 
primary importance. It can be easily modified to suit 
dynamic and stability analyses. . ; 

... | 
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CHAPTER 2 

P TRITE ELEMENT METHOD FORMULATION POR 
LARQ-E DEFLECTION ANALYSIS 

2.1 GENERAL 

The idealization of a continuum as an assemblage of 
individual pieces or elanents has often been used as a device 
for obtaining approximate solutions to problems which are 
insoluble in their original form (59). This discretization of 
the structure results in the replacement of the original 
governing ordinary or partial differential equations by a set 
of linear algebraic equations. With the general availability 
of high-speed electronic computers, this reduction of the 
problem is significant, since the solution of large systems 
of linear algebraic equations poses no difficulty. When two 
or three .i.dimensional elements are employed, the method is 
generally referred to as the ' finite— element ' method. 

The finite element method historically has been 
visualized by supposing the structure to be 'cut up' or 
divided into a number of sub- regions or elements, which are 
interconnected at a finite number of points, or ’nodes'. The 
force-deformation properties for an individual element, 
expressing relationships between nodal forces and the 
displacements are established usually by employing a 
variational principle. Once this relationship is estab- 
lished, the remainder of the analysis follows the usual 
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procedure of the stiffness method of structural analysis* 

Thus the d termination of the nodal force-deformation 
relations for the elemsit, or 'element stiffness matrix* 
which is the significant portion of the analysis. 

If the displacements of the structure are regarded as 
the fundamental unknown quantities, the principle of minimum 
total potential energy may be employed to obtain the element 
stiffness matrix. Other variational principles can be 
employed in the finite element method (60,20,61), however 
the potential energy principle is most often used (20,19) 
and is used here. 

In any numerical approach to a given problem 
convergence is very important. This problem of convergence 
in the case of the finite element method was investigated 

V X 

(20,60,61). Commendable work has been done by Oliveira (62) 
and Fried (63). Most of the early investigations of 
convergence simply involved successive refinements of grid 
size and comparisons with existing analytical solutions. 

While this procedure certainly does not prove convergence, 

it at least allows one to gain some confidence in the 

method in the absence of theoretical assurances of convergence. 

Melosh (20) formulated the finite element method in 
terms of the principle of minimum potential energy, and set 
forth some criteria for selecting displacement fields. He 
distinguished between three types of errors which can occurs 
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A) Idealization errors, such as those involved in 
replacing a shell structure by an assemblage of flat 
el em ent s . 

B) Discretization errors, such as those involved in assuming 
the form of the displacements within an element* 

C) Manipulation or round-off errors in the arithmetical 
operations. 

Variational principles can, at best, guarantee the 
convergence of the discretization error to zero as the grid 
size is refined. Thus, in order to relate the finite element 
method directly to variational methods, the idealization 
error should be minimized, if possible. The use of curved 
elements achieves this end in the present case. 

The admissibility criterion of the coordinate 
functions in the conventional Rayleigh-Ritz method is 
equivalent to the criterion of ’sufficient' continuity *of 
the assumed displacements, together with the imposition of 
the forced boundary conditions at a later stage in the 
analysis. Likewise, the usual ’completeness' criterion of 
the Rayleigh- Ritz method is analogous to the inclusion of 
the rigid body displacements and constant strain states in 
the assumed displacements. Melosh (20) and Bazely et al (6l) 
have enunciated these requirements. It must be clearly 


* 'Sufficient' continuity here means that u,v,w,w, x » and 
w,y are continuous over the entire structure. 
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understood that the above conditions guarantee monotonic 
convergence only to the true value of the potential energy. 

Uniform convergence to the true displacements and derivatives 
of displacements (i.e. stresses) does not necessarily follow. 
Negligence of terms representing a constant strain state has 
apparently resulted in convergence to incorrect results in 
applications to plate bending problems ( 63 ). 

Bazely et al (61) have presented the argument that 
the only necessary requirements for convergence to the true 
energy level are the inclusion of a complete rigid body 
motion and all constant strain states, although convergence 
is no longer monotonic in this case. The conclusion that 
these requirements are sufficient to guarantee convergence 
is reasonable since successive refinements of grid size lead 
to ever increasing satisfaction of continuity between 
elements, thus producing, in the limit, an admissible dis- 
placement field. In support of this argument, results of 
plate-bending analysis comparing 'conforming* and ’non- 
conforming' displacement assumptions were present (61 ) and 
showed, in almost all cases and for all grid sizes, that 
the so-called 'non^-conf orming' displacement field produced 
superior, although not monotonically converging, results. 
Before proceeding further some definitions are stated below 
for clear understanding. 
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Definition 1 ; Displacement models, displacement functions, 
displacement fields, or displacement patterns are the simple 
functions which are assumed to approximate the displacements 
for each element. 

D efinition 2 s 

w(X) = + a 2 X + X 2 + . . . + a n+1 X* 1 * 1 (2.1) 

Generalized coordinates or generalized displacement amplitudes 
are the a's, the coefficients of the polynomial in equation 
(2.1). The number of terms retained in the polynomial 
determines the shape of the displacement model, whereas the 
magnitudes of the generalized coordinates govern the ampli- 
tude. These amplitudes are called 'generalized' because they 
are not necessarily identified with the physical displacements 
of the element on a one-to-one basis; rather, they are linear 
combinations of some of the nodal displacements and perhaps 
of some of the derivatives of displacements at the nodes as 
well. The generalized coordinates represent the minimum 
number of parameters necessary to specify the polynomial 
amplitude . 

Definition 3 s 

u = a x + a 2 X + a 3 Y + a 4 X 2 + X I + Y 2 (2.2) 

The displacement model in equation (2.2) is expressed in 
terms of generalized coordinates, a , and is referred 
to as Generalized coordinate displacement model. 
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Defin ition 4 ? Conforming or compatible elements are those 
in which the displacement models are continuous within the 
elements and the displacements are compatible between 
adjacent elements. The compatibility condition implies that 
the adjacent elements must deform without causing openings, 
overlaps, or discontinuities between the elements. 

Definition 5 ; Complete elements are those in which the 
displacement models include both the rigid body displacements 
of the element and the constant strain states of the element. 

Definition 6 ; Geometric isotropy, spatial isotropy, or 
geometric invariance is that the displacement pattern should 
be independent of the orientation of the local coordinate 
system. 

Definition 7 ; The degrees of freedom of the element are the 
nodal displacements, rotations, and/or strains necessary to 
specify completely the deformation of the finite element. 

The degrees of freedom differ from the generalized 
coordinates in that each is specifically identified with a 
single nodal point and represents a displacement (or rotation 
or strain) having a clear physical interpretation. 

Definition 8 : The minimum number of degrees of freedom 
(or generalized coordinates) necessary for a given element 
is determined by the completeness requirements for convergence, 
the requirements of geometric isotropy and the necessity of 
an adequate representation of the terms in the potential 
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energy functional. Additional degrees of freedom "beyond the 
minimum number may be included for any element by adding 
secondary external nodes or specifying as degrees of freedom 
higher order derivatives of displacements at the primary 
nodes. The latter approach is preferred because it leads 
to a more compact numerical formulation in the assembly 
process. Elements with additional degrees of freedom are 
called higher order elements. 

Definition 9 : A natural coordinate system is a local system 
which permits the specification of a point within the element 
by a set of dimensionless numbers whose magnitudes never 
exceed ■unity. Moreover, these systems usually are arranged 
so that some of the natural coordinates have -unit magnitude 
at primary external nodal points. Not only does such a 
coordinate system generalize and simplify the foimulation 
but also facilitates numerical integration which is required 
to obtain element stiffness in linear case and tangential 
stiffness matrix in nonlinear (geometric) analysis. 

Definition 10 : An interpolation function, also known as a 
shape function, is a function which has unit value at one 
nodal point, and zero value at all other nodal points. 

Definition ~n : Isoparametric elements are those in which 
the geometry and displacements of the element are described 
in terms of the same parameters and are of the same order. 
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Defini tion 12 ; It is not necessary that the geometry and 
displacements of an element be expressed by the same order 
model. Subpar ametric elements are those elements for 
which the geometry is determined by a lower order model 
than the displacements. Elements for which the converse 
is true are called super par ametric elements. 

There are two types of boundary conditions, geometric 
boundary conditions or forced boundary conditions and natural 
or free boundary conditions. Geometric boundary conditions 
can be further categorized as being homogeneous or nonhomo- 
geneous and as normal or skewed. 

In practice, it may be true that neither inter-element 
continuity nor inclusion of a complete rigid body displace- 
ment is all important, as Haisler and Stricklin (64) have 
shown that omission of rigid body terms for shells of 
revolution does not affect convergence. This conclusion 
is borne out by results, not reported herein, of comparative 
studies using element stiffness matrices derived subsequently 
with and without inclusion of complete rigid body motions. 

A choice, arises between achieving inter-element continuity 
and including a complete rigid body displacement. One of 
the examples will show later on the results obtained 
choosing the latter. When the finite element method is 
regarded as a physical idealisation rather than a mathematical 
one, it is better to deal with a stiffness matrix which is 
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'equilibrated' that is, one which gives rise to self- 
equilibrating nodal forces due to nodal displacements, and 
this is achieved through the inclusion of a complete rigid 
body motion in the assumed displacements. In principle 
it is possible to include lagrangian multiplier terms to 
take into account the effect of discontinuities (64). 

However, ' this signif icantly increases the computational 
effort involved in obtaining a solution. It seems preferable 
instead to refine the grid size in order to obtain grea.ter 
accuracy. 


2.2 GEOMETRICAL RELATIONS 

C-f, ^ represent a set of orthogonal curvilinear 
coordinates for the middle surface and C the normal coordi- 
nate. The differential arc lengths along the parametric 
lines are denoted by dS. ( J =1,2) and the differential 

J 

surface area by dA. t 2 »n define unit vectors pointing 

in the ? 2 and r, directions . Referring to Figure 2.1, 

the geometrical relations are 

dS. = a. dS. 

3 0 3 



3 




n = t^ X 
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dA 


R. 




dS 1 dS 2 


oc. oc_. n r, ij 


(2.3) 


(i> j = 1,2) 


i 0 


in which r 


a . 


the position vector to a point on the middle 
surface, 

the surface metrics 


and 


R ij 


th e cur vat ures. 


For the cartesian case, becomes X and £ 2 becomes 
Y. The surface, Z = Z(X, Y), is said to he shallow when 


Z,/« 1 ; Z,y 2<< l 

Specializing (66) for the shallow case leads to 


( 2 . 2 ) 


R 


a l ~ a 2 ~ 1 
\ * X X + Z 'X *z 

= ly + Z,y 

n = “Z , y iy ” Z , y iy + i z 

T 

-Z,. 


(2.4) 


11 


XX 


... 3 : .... gs —Z. 

R " »yy 

2 ^ ■ . ****** 




12 


”Z , 


XX 
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2.3 deformation-displacement relations 

Referring to the Figure 2.2, the strain matrix [e] 
at- (H^, S 2 » c) is expressed as 


[ e ] ={ v e 2 , Ei2 } 

2 9 ^ 


( 2 . 5 ) 


in which 


extensional strains ( ^ » directions); 


and e^ 2 = the shear strain. Considering the shell to he 
thin and neglecting transverse shear deformation, the 
strain varies linearly over the thickness and can he 


expressed as 


[e] = [e] +C [k] 


( 2 . 6 ) 


in which [e] contains the stretching deformations and [k] 
contains the curvature changes for the middle surface. 

The expressions for [e] and [k] used in this study 


[e] = [e] o + Cel 


U 1,X “ Z, XX W 

Ug *y 2 f yy "W 

U 1,Y + U 2,X " 2 Z, XI w 


[e] r = ijw 


(2.7) 


| 2w, x w» T 
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[k] 


w. 


w , 


XX 


'YY 

2w, xr 


Pi = -w, x 
P 2 " W, Y 

These relations are restricted to small rotations, i.e. 

1 . 

2.4 FORCE-DEFORMATION RELATIONS, STRAIN ENERGY DENSITY 

Stress matrix [o] corresponding to [e] is expressed 

as 

[a] = {a lf a 2 , a^y (2.8) 

and express the stress-strain relations as 

[d] = [D] <[e] - Ce]°) ( 2 - 8 a) 

i» *1 0 

in which LeJ contains the initial strains. It is assumed 
that initial strains also v-ry linearly over the thickness, 

' t ' , and write 

[c)° = [e)° + C[k]° 

in which 

[e] ° = | I [e]° d? 

[k]° = ^ /?[a]° d? 
t 3 

The expressions for the stress resultants and stress couples 
follow from the definition, expressed as 


(2.9) 


(2.10a) 
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W = J O] dc = { N r n 2 , n 12 > 

[M] = fc O] dc = { M x , M 2 , M^} 

and are written as 

[K] = [H]° + [B] s ([e] t + [e] r ) 

[M] = [M]° + [B] b [k] 


in which 

[H]° = -Mg [ej° 


and 

[M]° = -M b [k]° 


[D]g = = Stretching rigidity matrix 

t 3 _ 

[D]^ = “igl-D] = Bendin g rigidity matrix 


(2.10h) 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


V represents the strain energy per -unit area of 
the middle surface. By definition, the first variation 
is given hy 


6V =/ [of 6[e] d$= [N] T 6[e] + [M] T 6[e] (2.14a) 

The second variation is obtained hy operating on 6V, 
expressed as 

6 2 V = 6(6V) 

= [N] T 6 2 [e] + 6[N] T 6[e] + [Mj T 6 2 [k] + 6 [m] T 6[k] 

(2.14h) 

U.T. KANPUR I 
central libr ary 5 

ACC. No. A f 189.1 | 
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As [e] and [k] are linear functions 
6 a [e] * = 6 2 [k] = [0] 

6 2 [e] = 6 2 [e] (2-14o) 

Substituting for [N] and [M] using equations (2.11) and 
noting that [N]° and [M]° are constant, 

6V= (6[e] 4 + 6[e] r ) T [[H]° + [D] g ([e] ^ + [e] r )] 

+ (6[k]) T ([M]° + [D] b [k]) (2.15) 

6 2 V= (6 2 [e] r ) T [1ST] + (6[e]j + 6[eJ r ) T [D] g (6[e] 1+ 6[e] r > 
+ («[k]) T [D] b 6[k] (2.16) 


2.5 ELEMENT STIFFNESS AND NODAL FORCE MATRICES 

[u] and [p] define the middle surface translation 
and external force intensity matrices referred to the initial 
local reference frame as shown in Figure 2.2, expressed as 


[u] = { U-^ n , w> 

(2.17) 

Cp] = { Pq» P 2* P n } 

The governing equilibrium equations are obtained by 
applying the principle of virtual, displacements which 


requires 

/ 6V d A = / [p] T 6[uJ d A 


(2.18) 


to be satisfied for arbitrary permissible 6[u]. Also, the 
consistent linearized incremental equations are determined 
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By taking variation of equation (2.18), expressed as 

/ 6 2 V a A = / 6[p] T 6[u] d A (2.19) 

In this section, the expansion of equation (2.18) is 
considered. Equation (2.19) is utilized later on in 
connection with the Newton-Raphson method. 

The expansion of [u] over element n can he expressed 
as 


[и] = [A] [a] n (2.20) 

in which [a] contains the generalized coordinates for 
element n and [A] contains prescribed functions of the 
local coordinates. 

Substituting for [u] in equation (2.6) leads to 
Ce]„= [B]„. 

( 2 . 21 ) 

XI 

[к] = [B] b [a], 


l 


SI 


[e] r = [B] sr [«] n 


'n 


The variations are 

«[e] t - [B] Bi «[«] B 

«[e] r = [B]* r 6[a] n (2.22) 

6[k] = [B] b 6[«] n 

It is significant to note that [B]g r and [B]g r depend 
on [a] whereas [B] o0 and [B], are independent of [a] . 
Substituting equations (2.21) and (2.22) into equation (2.15): 
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and integrating over the element* the result may he written 
as 

/ 6V dA = 6[«]* [ [E] n C«] n + [p]° ] = 
element n 

(2,23) 

in which 

[K] = / ^ B ^Sr^ ^S&^^S + ^Sr^ 

element n *• 

+ [B] T [D] h Cb] ] dA (2.24a) 

b u u 

and 

[S]° = I t([B] +[B]* r ) T [n]° + [B]^[M] 0 ]dA 

element n 

(2.24b) 

The term involving the external loading is written as 

/ [p] T 6[u] dA = 6[a] n ^] ext>n (2.25) 

element n 

in which 

[B] ext,n = / [A]T ^ dA (2.26) 

el men t n 

2.6 LINEAR STATIC CASE 

Eor the case of the shell undergoing time- 
independent deformations, Hamilton's principle becomes the 
principle of minimum potential energy. If the deformations 
are assumed small -linear static case - [B]^ and [B]^ 
matrices are zero. Thus, the element stiffness and force 
matrices are. 
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element n 

/ 

[£]° = / ([B]g £ [N]° + [B]^' [M]° dA (2.27) 

element n 

[P] ext,n = 1 [A]I [p] dA 

element n 

The general system of equations is, 

[rig, [u] = [P] (2.28) 

The displacement boundary conditions can be applied to 
the above and the resulting system of linear simultaneous 
equations can be solved for the nodal displacement vector 
[u] . Once [u] is known, the element nodal displacements 
[u] n for each element can be obtained. Thus using equations 
already shown, strains and the corresponding forces and 
moments can be obtained. 

2.7 GEOMETRICALLY E0N1IKEAR STATIC CASE 

This case corresponds to one of time- independent 
loading on the shell but with large out of plane rotations. 
In this case, the nonlinear terms inte^have to be taken into 
account. Thus, the terms, [B]g and [®]g r » are & ue 
geometrical nonlinearity. The expressions for [e] r and 
6[e] r are used to generate the matrices [®]g r an( ^ [®]g r * 
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. <W ’Z )! 

[e] r - [B] gr [«] n - + ^ . J. 

2W f Wfy | 

w, 2 6w, z 

«[®] r = CB]* [a] - j w-y 6w ’T 




w,j 6 w, y + w, Y 6w, 


Z 


J 


(2.29) 


The form of equation (2.29a) suggests working with 
{w » yr> w, y } and 


!>]*= 


W 'X 0 


o 


w 


*Y 


w f *y" ^ j 


(2.29b) 


Then [e] r = t [«]‘ 


w 


’Z 

w, T j 


. j 6w * X j 

6 [e]= O] f 

1 5w »yJ 


Continuing 


(2.29c) 


W, Z 


W 


= [A] r C«] n 


'Y 


(2 .3*0) 


It follows that 


[E] Sr - 

[®] g r = [“]'[A] r 


(2.31) 
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[ w] * depends on [a] and the local coordinates. It is 

convenient to segregate the nonlinear terms in the expression 

0 


• for [K] n and [P] n . Therefore, 
[K] n - + ^r,n 


[5]° = [P]? + [p] u 

Lrj n J l,n J r,n 


«-,0 


( 2 . 32 ) 


in which [K] n - » [Pj? _ are the linear terms, expressed 

Jo | XI Jo f XI 


as 


[K] 4n = / ([B]^[D] g [B] sl+ [B]“ [D]„ [B] t )4A 

’ element n ^ ^ 

[F]°, n = / <[B 4 M° + [B]* M°) dA 

element n 

and [K] r>n and [P]° >n depend on [a] n> expressed as 

[g h,n= , ! + tD] S [B C + ([B] sJ [:D] S [B 4 )T 

’ el ement n x ' 

+ «tB]| r ) T [D] g [B]* r ]dA (2.34) 

[P]° >n = ! [K]°]dA 

element n 

Explicit forms for the integrals can be obtained for the 
nonlinear matrices but the computations are difficult. 
Pellippa (40) and Gallagher et al (84) assume [ w] is 
constant over the element. However, it is preferred to 
use the exact expression for [w]* and evaluate the 
integrals numerically. 
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It remains to substitute for the element parameters, 
[a] n , in terms of the nodal displacements. Five degrees of 
freedom per node are considered, namely, middle surface 
translations u^, u^ and w, and normal rotations and p^. 
The number of nodal displacement quantities depend on the 
expansions for polynomials used. For example, Schmit (86) 
worked with 12 per node. The above applied both for 
nonlinear flat plate analysis and for shallow shell analysis 

[U]. denotes the displacement matrix for node j, expressed 

J 

as 


[u] = { , u 0 , w, p , p } (2.35) 

3 1 d id nQde j 

[U]_ represents the matrix containing the nodal displace- 
ment matrices for element listed in a prescribed order. 


[U]_ ={ [U] n , [U] 


E,n 


J n, 


frh } 


J n 


Sn 


( 2 . 36 ) 


in which ^ = the total number of nodes and n^ (i=l, 2, . . . , 
are the node numbers. Figures 2.3 shows the notation 
specialized for a rectangular element having nodes at the 
corner points. 

Using equation (2.20), 


[U]_ = [C] [a] 

L J E,n J n J n 

[a] n = [Cj; 1 [0] Ein 

Substituting for [a n ] in equations 


(2.37) 

(2.23) and (2.25), 




FIG 2*3 ELE MENT NOTATION: g^ft|SIAI>l 



60 


element n 


/ 67 dA = 6[U]£ jn ([K] n [0] Bjn + [P]°) - «K, n [P] intjn 


/ [p] T 4A - «[*]£,„ [P] extjn 
element n 


,T 


( 2 . 38 ) 


in which 

[ K ] n = [c ]; 1 - 1 [i] a cc ]- 1 

[ p ]n = tP]° (2-59) 

Wert.n " C °]" 1,T t 5 hxt,n 

interpreting* [K] n and [P ] as element stiffness and nodal 
force matrices. Finally, substituting for [K] and [P]: 
using equation (2.32) 


and 


# 




+ [k] 


r ,n 


[P] 


0 

£ ,n 


+ [P] 


0 

r ,n 


( 2 . 40 ) 


are obtained, in which the terms in [K] and [P]^ are 
determined by applying equation (2.37) to equations (2.31) 
and (2.32). 

2.7.1 System Equations 

[u], defined as the nodal displacement matrix, is 
expressed as 

[u] ={ [n] 1 , [u] 2 [U] } (2-41) 

in which n^ represents the total number of nodes and the 
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summat ion of equation (2.38) over the elements is expressed 
as 


n 


E 

l 

n=l element n 


/ 6V dA = f ^E.n [P] lnt = S[u] T [p] 

**T 


n=l 


n 


E n E 

I / [p] T 6[u]dA = I 6Wj n [P] ext n = 

n=l element n n=l h,n ext,n ext 

Applying the Principle of Virtual Displacements to equation 
(2.42) results in equilibrium equations 

[ *] = [p] ext - [p] int = [+] (I'll) = [0] ' (2.45) 

In this study, the external loading is considered to be 
independent of the displacements, i.e. take [p] ex ^ = constant. 
Term is a nonlinear function of [u] when the 

rotational terms are included. 

2.7.2 Newton-Raphson Iteration 

In the following it is assumed that the nodal dis- 
placement boundary conditions have been introduced. 

[u] (i) represents the ith estimate for the solution 
due to a particular loading. 


(2.42) 
ilr 


(i) 


[#] lw= w w = [ * ] 


(i) 

Cp]lnt l[u]=[ u ] (i) ” p]lnt 


(2.44) 


The recurrence relation for Newt on-Rap h son iteration is 
-6[>] (i) = 0] (i) (2.45a) 
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in which is the first order change in [ip] due to 


6[u]. Equation (2.45a) "becomes 


= ^ext 




(2.45b) 


6[p] ini can be expressed as 


6[p] ^nt = Wf 5 6[u]<i> 


(2.45c) 


in which [k]!"^ can be interpreted as the system tangent 
T (i) 

stiffness matrix evaluated at [u] v . Finally the 
recurrence relation takes the form 


M 1} { M (i) - [pw - 


and [u] (i+l} = [u] (i) + 6[u] 


( 2 . 46 ) 


A convenient measure of [u] is the Euclidean norm, E, 
defined by 


N = ([u] [u]) T 


(2.47) 


The convergence criteria is 


r(i+l) _ w (i) 


( 2 . 48 ) 


in which A = the desired accuracy. 


2.7*5 Generation of Tangent Stiffness Matrix 


By definition 


-f / 


6[u] x [k], 6[u] = l f 6 2 V dA 

n=l element n 


(2.49a) 
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Prom equations (2.16) and (2.22) 

/ 6 2 V dl = / [(6 2 [e] ) T [N]dA 

element n element n r 

+ 6Ca] “ { element ‘ (CB] B + +&V 

+ l 3]J [D] t [B] b ]dA} 6[a] n (2-493) 


To proceed further, the first term 


(*' S [e] r ) T [N] = + N 2 (6 w, y ) 2 + 2N 12 «5w, x 6w, y 

(2.49c) 

must he expressed in a quadratic form in 6[a] n - By 


defining [Z]* as 
[N]* = 


»1 N 12 


N 12 N 2 


(6 2 [e] r ) i [H] = 6[a]£ ([i] r W*[A] r ) S[a] n 


■Ir.i -r ,T 
*n 

can he written. Then 


T 


(2.49d) 


(2.49e) 


/ 6 2 ydA = 6[a]hE], r 6[a] 
element n n x,n 


in which 


[Zl . = [Z] „ + [f] 

L J t,n L J £,n J g,n 


n = 6 tn]^ n Cz]^ n 6[u] E#n 
(2.50) 


Z 




(2.51) 



The linear term is given by equation (2.33)* The nonlinear 
term, which is usually referred to as the geometric stiffness 

matrix, has the form 


[El = I [ULVi'IaL + W'.WJjL 

g ' n element n r r b b bT 

* ([B]^ CB] S [B]^.) T + ([B]* r ) T [B] s [B]^]dA 

(2.52) 

It may he noted that [K]„ involves the same terms as 
[K]„ the only additional term is ,{A]^ [U]*[a] * 

X f xl X X 


2*7.4 linearized Incremental Procedure 


[u] denotes the solution for a partucular loading. 
An incremental loading, A[p] ^ is applied. The incremer 
equilibrium follows from equation (2.43)* expressed as 


A ext A ^-^int 


To determine A [p] . 


(2.53a) 


int' 


«[u] T A[p], + = f / A(«V)dl 

n=l element n 


n- 


E 


(2.53b) 


Taking a(6V) - 6 2 V results in linearized incremental 

equations 

i[ P]int “ 6[p] int = [K] t a[u] (2 - 53c) 

The incremental displacement, A [u] , due to 4[p] ex ^. is 


obtained from 
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in which [K]^ is evaluated at [u], i.e., at the equilibrium 
position prior to the application of the load increment. 

The procedure corresponds to truncating Newton- 
Raphscn after one cycle. As the solution is approximate, 
the only resort is to repeat the computation using smaller 
load increments in order to assess the accuracy. 
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CHAPTER 3 

SOLUTIONS POR SMALL DEFLECTION 


3.1 GENERAL 

The formulation outlined in the previous chapter is 
valid for an arbitrary element, i.e., arbitrary shape and 
displacement expansions. A simple shallow rectangular 
element as shown in Eigure 2.1 is employed. The displacement 
expansions are 


u 

v 

w 


ol, + a_X + a„Y + a.XY 
12 3 4 

(3.1) 

a + a,X + a Y + aXY 
3678 

a 9 + a io X + a ll Y + a 12 X2 + a l3 XI + a 14 Y2 

+ a 15 X 3 + a^X^Y + a-^YY 2 + a i8 y3 + a i9 X3y 

+ a 20 ^ 3 ‘ 


The element is unconforming (it does not satisfy 
normal slope compatibility) and complete rigid body motion 
is not taken into account. This point has been discussed 
in Chapter 5. In spite of non- compatibility and ommission 
of complete rigid body motion, it has been shown (82,64) 
that results obtained are still accurate. Hence the above 
simple displacement expansions are chosen which still give 
accurate results. Linear stiffness matrix for shell element 
is generated using numerical integration procedure as outlined 
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in Appendix A, as stiffness matrix given by Brebbia ( 67 ) 
failed to give correct results. Further the stiffness 
matrix obtained by numerical integration technique is checked 
for its accuracy by using the technique su^ested by Kabaila 
as given in Appendix B by writing the necessary computer 
program. 

The above developed stiffness matrix has been used 
in the program of the finite element method to analyse 
hyperbolic paraboloid shells bounded by straight edges with 
different boundary conditions. Since displacement formula- 
tion is used in the finite element theory by applying the 
Principle of Minimum Potential Energy, only forced boundary 
conditions are specified in the analysis. The natural 
boundary conditions, together with the equilibrium equations 
are satisfied approximately as part of the minimization 
procedure. Since forced boundary conditions are specified 
in the analysis, it might, therefore, be expected that 
convergence would be superior for those problems in which 
all or most of the boundary conditions are of the forced 
type. Also, the extent to which the natural boundary 
conditions are satisfied for a given problem can be used 
as a means of assessing if convergence is satisfactory in a 
practical sense. 

A variety of problems of hyperboloid 
paraboloid shells bounded by characteristics for which 
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analytical or other numerical solutions are available, are 
analysed. These problems are intended to check the computer 
program developed independently and to substantiate the 
finite element method as applied to shallow shells and' to 
enable one to evaluate in a limited sense the convergence 
characteristics of the element stiffness matrix generated. 

3.2 SIMPLY SUPPORTED HYPERBOLIC PARABOLOID SHELL BOUNDED 
BY CHARACTERISTICS UNDER UNIFORM NORMAL LOAD: 

COMPARISON WITH SERIES SOLUTION (7) 

Figure 3*1 shows a square hyperbolic paraboloid 
simply supported on all edges and subjected to a uniform 
normal load of -50 lbs/sq.ft. 

Shell Data: 

a = 180 in. E = 3x10^ lbs/sq.in. 

b = 180 in. v = 0.16 

c = 36 In. 

h = 2.5 in . 

A grid of (8x8) over the entire shell is used in the 
analysis. 

Displacement boundary conditions for an edge X — constant 
are : 

v = 0 
w = 0 

e 2 - o 




(b) Variation of Mx along x axis 



(cl Deflection along x axis Mes n size 8x8 

Data E = 3x10 lb/ in 2 
a =180* 

p >50 in/ ft 2 
c=36 ■ N ^ 

h*2-S" 
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Similar boundary conditions are used for an edge 
Y e= constant. The normal deflection w, the membrane shear 
Nyy and the bending moment My across the section Y = 0 is 
plotted along X-axis from Z = 0 to I = a. 

All quantities plotted agree with the series 
solution. The membrane theory predicts a constant Nyy 
of 156.25 lbs/in. throughout the shell, while the bending 
theory yields a value of Nyy = 160 Ibs/in. at the centre 
of the shell. The bending moment My is fairly large, 
producing bending stresses which are about half as much as 
the membrane stresses, and dies out slowly toward the 
centre of the shell. 

3.3 CLAMPED HYPERBOLIC PARABOLOID B GUIDE D BY CHARACTERISTICS 
UNDER UNIFORM NORMAL LOAD 

3.3.1 Figure 3*2 shows a square hyperbolic paraboloid 

clamped on all edges and subjected to a uniform normal 
load of -1.0 Ibs/sq.in. This shell has been analysed by 
Chetty and Tottenham (7) using a variational method. 

Shell Data: 

a = 6.46 in. E = 5x10 lbs/sq.in. 

b = 6.46 in. v = O .39 

c = 1.304 in. 
h = 0.25 in. 

A grid of (8x8) over the entire shell is used in the 
analysis. 



04 * ( uO M 1° 
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Displacement "boundary conditions for an edge X = constant 
are: u = 0 

v = 0 
w = 0 

"i- 0 

P 2 = 0 

Similar boundary conditions are used for an edge 
Y = constant. 

The figure shows the comparison of the finite element 

analysis with those given in (7). This agreement is quite 

good considering the size of the grid used in the finite 

element analysis. This is in agreement with the fact that 

there would be improvement of accuracy when the boundary 

conditions as in this case are all of the forced type. It is 

interesting to note that the agreement is better for the 

membrane shear force N than for the bending moment M Y . 

XY x 

The membrane shear N n exceeds the membrane value of 

16.0 lbs/sq. in. by only a slight amount at the centre of the 

shell , but the bending moment propagates into the central 

region of the shell, producing fairly large bending stresses. 

The maximum deflection at the centre of the shell is 
-3 

8.732x10 in. The maximum deflection at the centre of the 
shell with knife edge supports (26) is 9.18x10“^ in. The 
displacement boundary conditions at I = constant are : 


u = 0 
w = 0 
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Similar boundary conditions are used for an edge 
Y = Constant. It is seen that the central deflection of a 
knife edge supported shell is about 1.14 times that of the 
shell with all edges, clamped whereas inplates it is about 
3.3 (85). 

3.3*2 Figure 3*3 shows a square hyperbolic paraboloid 
clamped on all edges and subjected to a uniform normal load 
of -0.1 Kg/sq.cm. This shall has been analysed by Connor 
and Brebbia (24)- 

Shell Data: 

a = 50 cm 
b = 50 cm 
c = 10 cm 
h = 0.8 cm 

A grid size of (8x8) over the entire shell is used in the 
analysis. 

Displacement boundary conditions for an edge X — constant 
are : 


u = 0 
v = 0 
w = 0 



E = 28,500 Eg/sq.cm 
v = 0.4 


0 




Finite element. { 8 x S }■ 
(Exactly same; as brebbia } 
Finite difference HS x 16 } 
Kantorovich solution 




JORMAL DEFLECTION AT CENTER LINE CL Ah' FED HYPAR 

(fsfOER: U niform normal load 



/ / 


Mesh size 8x8 


a? by cm 


b=50cm 
E- 28,50 


cm 


F^<)lkgcTh 2 
c = 10 cm 
h = 0*8 cm 
v = 0*4 
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Similar boundary conditions are used for an edge 
Y - constant. 

The Figure 3 . 4 shows a comparison of the finite 
element results for the normal displacement at the centre 
with the values obtained using the Kant orovi ch-G-al erkin 
procedure (7) and second order finite differences (67). 

The variation of normal displacement at the centre- 
line (Y = 0) with mesh size is plotted in Figure 3*3* Close 
agreement is obtained with a finite element mesh size 
(8x8) equal to twice the finite difference grid spacing 
(16x16), the percentage difference is less than 4 percent 
for the central deflection. By finer mesh size, the 
finite elemaat solutions may be made to agree within 
1 percent at all points. 

The comparison of results shows that one can obtain 
good results with the finite element method using a 
reasonably coarse mesh and fairly simple displacement 
expansions which are not compatible and which do not 
include all the rigid body displacement terms. 
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CHAPTER 4 

SOLUTION FOR LARGE DEFLECTION OF CLAMPED HYPERBOLIC 
PARABOLOID SHELL BOUNDED BY STRAIGHT EDGES 

4.1 APPLICATION TO HYPERBOLIC PARABOLOID BOUNDED BY 1*3 
. CHARACTERISTICS 

In order to demonstrate the capability of the finite 
element formulation presented earlier, a clamped hypar shell 
bounded by its characteristic lines is chosen. The Newton- 
Raphson method has been employed to trace the load-def lection 
curve. Since the computer memory is limited, a 4x4 mesh 
size is taken. The tangent stiffness matrix is evaluated 
by numerical integration method using Gauss-Legend re 
quadrature. The load-deflection curve is drawn for the 
centre of the shell. The nondime ns i on al parameters used 
are 

P = a(l2(l- 2 )) 1 / 4 /(4tR) 1/2 

w = w/t ( 4 .1) 

Q = P N R 2 /Et 2 

v = 0.4 

where a is the dimension of the square planform along 
X-axis, R is defined by the surface equation of the shell 
Z = XY/R. 
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For the shell under consideration 


R = = 250 


0.8 


v = 


R _ 

t 

H = 


0.4 

512.5 

5-15 

Pjj R 2 /st : 


= 5.247 E 


N 


4*2 DISCUSS ION OF THE RESULT 

For the buckling of the shell, experimental 
results as shown in Figure 4.1 give a value of 0.558 Kg/ cm 2 
(6*7). Buckling is considered to occur when the behaviour 
changes abruptly, i.e., when there is a sudden change in 
the load-de flection behaviour. Figure 4.1 Indicates 
buckling load as 0.48 Kg/cm 2 . The only rigorous mathe- 
matical equation for predicting the buckling load of an 
ideal hyperbolic paraboloid, simply supported on the 
perimeter, was derived by Reissner (5 )» namely 


2Et 2 h 2 


N cr 


5(l_v 2 ) 1 / 2 a 2 b 2 


(4.2) 


where = the magnitude of the uniformly distributed 

load that causes buckling. Equation (4.2) is most 
significant because it establishes the functional 
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relationship between the buckling load and the pertinent 
shell and material variables. The form of the half-wave 
length of buckle, A , is given by 

T) 


A 

P 


2 (ab) 1 / 2 ( | ) 1/2 

2 4 0(l-v 2 )} 1/2 


(4.5) 


and is physically interpreted by Reissner to represent a 
'wrinkling' of shell in the direction of the compression 
diagonal. Reissner 's analysis based on simply-supported 
boundary conditions gives (lower critical load. The ex- 
perimental, finite- element and theoretical (Reissner* s 
linearized stability analysis) values of buckling load 
are compared in Table 4.1. 


Table 4.1 


Experimental (6 7) 

Finite element 
(4x4) 

2 

Reissner (5) 

1 

5 

0.558 

0.480 

0.580 


The experimental result gave lower value due to imper- 
fections and the probable lack of complete translational 
fixity at the boundary. The computation has been done on 
IBM 575/155 Computer at I.I.T. Madras. 
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4.3 EXPERIMENTAL CONCLUSIONS OP LEET (68) 

Some of the experimental conclusions of Leet (68) on 
the study of stability in the hyperbolic paraboloid bounded 
by characteristics are mentioned below. 

A) Buckling of a fixed edge hyperbolic paraboloid 

is characterized by the entire surface of the shall forming 
ripples in the direction of the compression diagonal. 

B) At the critical load, the hyperbolic paraboloid 
surface undergoes a sudden bend or knee in the load- 
deflection curve. 

C) If sufficient edge member stiffness is present to 
prevent buckling of edge members before the shell buckles, 
it seems that the buckling load that a shell can develop 
may be related to the edge member axial stiffness withfna 
limited range. 

D) The hyperbolic paraboloid with fixed edge boundaries 
and loaded by air pressure, buckles into half-waves whose 
half-length is closely predicted by equation (4*3)* 

E) Imperfections reduce the magnitude of the initial 
buckling load. In the air pressure tests, the buckling 
load varied from approximately 51 percent to 71 percent 
of the load given by equation (4*2). The larger the 
imperfections, the lower the initial buckling load. 
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F) After initial buckling wave formation has occurred, 
the shell develops a tension field mode to carry additional 
load . 

G) The ultimate load of the shell and edge member 
system is shown to be dependent on the edge member moment 


of inertia, 1 -^ , P rov ' ided sufficient area is present to 
prevent fielding of the material and no local instability 


occurs. 


H) The ultimate load of a hyperbolic paraboloid is 
dependent on the manner of support of the edge beams. The 
shell whose edge members are fixed at one end and simply 
supported at the other end carried more load than the 
shell whose edge members are fixed at one end only. 
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CHAPTER 5 

CONCLUSIONS AND DISCUSSION 

5.1 CONCLUSIONS 

The finite element as herein applied to shallow 
shells, using elements of the same shape as the shell 
middle surface rather than a series of flat elements has 
been shown to give results in good agreement with exact 
solutions . The most advantageous feature of the method 
is its versatility and the ease with which, for example, 
various boundary conditions (including elastic edge beams) 
can be included. 

It is found that inclusion of complete rigid body 
motion in the displacement function makes no change in the 
results obtained in the case of hyperbolic paraboloid shells 
bounded by characteristics with clamped edges from those 
obtained without inclusion of complete rigid body motion 
(4). 

It is also reported (82) that Birkhoff-Garabedian 
interpolation formula which ensures compatibility i.e., 
the continuity of derivatives of the normal displacement 
across inter-element boundaries has not been essential 
as the results agreed with those obtained violating 
compat ibili ty . 
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Some of the peculiarities of the finite element 
results, notably the under-estimation of the membrane 
forces near the boundary, may be traced to the fact that 
the natural boundary conditions are not exactly satisfied 
for a finite grid size. In a problem for which the exact 
solution is not known, the degree of satisfaction of the 
natural boundary conditions may be used as a practical 
method of assessing whether or not a fine enough grid 
has been used in the analysis. 

5.2 SCOPE FOR FURTHER WORE 

1) The improvement of the satisfaction of the natural 
boundary condition involves displacement shapes which are 
capable of satisfying, for example, the membrane natural 
boundary conditions. This, however, would result in an 
increased nuinber of generalized coordinates per element 
and thus for the same computational effort, not as fine a 
mesh could be used. Therefore, whether or not improved 
accuracy should be sought by using more redefined displace- 
ment assumptions, or by using relatively crude assumptions 
together with fine grids is answerable only by actual 
numerical comparisons. 

Alternative methods of satisfying the natural 
boundary conditions are the use of Lagrangian multipliers 
(64) or a reformulation as a mixed variational problem in 
which some of the displacements- and some of the stresses 
are . assumed in terms of generalized coordinates (58, 7^). 
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2) Several extensions of the scope of this investigation 
suggest themselves. Temperature effects, non-uniform thick- 
ness, anistropy, and flexible column supports could all be 
immediately included. 

Other extensions, such as to creep, vibrations, 
buckling (both pre- and post-) and plasticity could also 
be made. 

3) Skewed and tilted hyperbolic paraboloid shells 
should also be brought within the purview of finite element 
method. 

4) It is necessary to develop a deep shell element 
using discrete Kirchhoff assumptions and shear strain 
displacement relations. Noteworthy start has been made 

by Gurbachan Dhatt and Jean-Louis Batoz (87) by developing 
a deep triangular shell element. Deep shells may necessi- 
tate the consideration of wind loads in addition to 
gravity loads. 

5) Thick shells are still a challenge. The straight- 
forward use of the three-dimensional concept offers certain 
difficulties. In the first place the retention of three 
degrees of freedom at each node leads to large stiffness 
coefficients for relative displacements along an edge 
corresponding to the shell thickness. This presents 
numerical problems and may lead to ill-conditioned 
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equations when shell thicknesses become small compared 
with other dimensions in the element. The second factor 
is that of economy. The use of several nodes across the 
thickness ignores the well-known fact that even for thick 
shells the 'normals’ to the middle surface remain practi- 
cally straight after deformation. Thus an unnecessarily 
high number of degrees of freedom has to be carried, 
involving penalties of computer time. Shear deformations 
are an important feature in the thick shell situations. 
Special isoparametric curved hexahedron three-dimensional 
elanent has to be fully developed and exploited. They are 
already proved to be better than tetrahedron elements. 
Obviously numerical integration methods come in handy. 

For axi-symmetric shells the formulation is obviously 
simplified as the mid-surface is defined by only two 
coordinates p and a considerable saving in computer 
effort is obtained. Spherical domes, edge loaded cylinder, 
cooling towers and curved dams require the thick shell 
analysis. 
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APPENDIX A 



Figure A.l: Normalized Coordinates for a 
Rectangle. 

u = X + Y + XY (A.l) 


Substituting the nodal coordinates in the above 



1 

-a 

-b 

ab 


1 

a 

-b 

-ab e 

{u> e = 




{a} 


1 

a 

b 

ab 


1 

-a 

b 

-ab 

= [a 

] {a> e 



(A.2) 

u = [l X Y X Y] 

[A]" 1 

{u} e 

(A-3) 


Simplyfying the R.H.S., 
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u . [ 


4 v*- f)d- 1) + fx 1 -# h 1+ i )(1+ 1 } 


a' v V 4 
Zw , 1 '-* - ’ e 


|(1 - |)(1 -+f)] {u> 


(A. 4) 


Let 5 


X 

a 


Y 

* = b 


u „ [Id-jjd-n) |(i + e)d-n) |d+c)d+n) fd-s) d+'i)]t » 1 

(A.5) 


Similarly 


v . [i (1 ^ )(!-„) i ( ne )d-n) j(i+?)d+n) id-e)d+r,)]{u}' 


(A. 6 ) 


w = A- 


1 +A 2 X + A^Y + A 4 X 2 + A 5 XT + A g Y 2 + A ? X 3 + AgX 2 Y 
9 — ' "iu 11 


+ A„H 2 + Ai 0 T 3 + An X 3 Y + A ^2 XI' 


1 3w 


^1 a 3 X 

= - l/a[A 2 +2A 4 X+A 5 Y+3A 7 X 2 +2A 8 XT+A 9 Y 2 + 3A 1IL X 2 Y + A^Y 3 ] 

(A. 7) 


6 = _ 1 JJi 

" 2 b 3 Y 


= - i[A 3 +A 5 X+2A 6 Y+A 8 X 2 +2A 9 XI + 3A 1() Y a + A^X 3 +3A 12 Efa ] 

Substituting the nodal coordinates of the element xn the 
above, 


U) e = [C] (A) 


(A. a) 
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or 

{ A } = [c ]"* 1 { w } 0 (A. 9 ) 

. * . w=[l? rj f £ p p 2 ^ n 2 r) 3 ^3 r] g p 3 ] [c]~ 1 {w} e 

or 

W = [% N 2 n 3 n 4 . n 5 n 6 f 7 n 8 S o N 1:L jsr 12 ] m e (a. 11 ) 

where to are polynomials in ,r , p . 

For a two-dimensional strain-displacement relationship, 
d 13- 

we need to know -g— ^ etc. These are not directly available and 
hence a relationship between X,Y and Z » P must be established. 
By partial differentiation 


3 u ^ 


: jl! 

3 Y 


' 3u N 

dZ 


3? 

3? 


9X 


L ~ 

f 




I . 

3 u ! 


J_X 

, 9 Y 


3 u | 

3 T) 

> 

> 

__ 3 P 

3 p _ 


3 Y ) 


a u ^ 
[J] j 3 X 
1j_u 5 

1st j 


where [j] is aJacobian matrix 


(A. 12) 


3 u ^ 


j 3 u~j 

3 X i 


■i 35 

J 


= [j]' 1 ! 

1 3u 


| 3 u j 

3 Y 

k t 


U p j 


(A.13) 


The relationship between the normal and global cartesian 
coordinates is given by: 
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x = + i-d+5 ){l-n)X g + i(l+E)(l+l)X 3 

+ l(i-«(i+n)x 4 

T = |(l-C)(l-tl)T 1 + j(l+5 )(l-l)T 2 + |(l+?)(l+ri)Tj 

+ A(l-?)(l+r,)7 ( . 

The Jacobian matrix for the above becomes 

f a/2 0 1 



[J] 


-1 


2/a 0 

0 2/b 


or 


’ 3 jh 


( 

2 

•\ 

3 X 



a 

a? 

\ 


< 


- > 

3u 

| 


2 


3Y 



Id 

^ j 


{ \ 
3 


2 _L_ 

3X 


a H 

JL 


, v 

2 _3_ 

3Y 


Id ^ 



J 


(A.14) 


(A. 14) 


(A«15) 


(A. 16) 


(A.17) 


(A. 18) 
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' 



B 2 w 

S x s 


a 2 

35 2 

3 2 W 

► = 1 

4 

djw [ 

1 3T 2 

r 

b 2 

f ■ ■ 

. ^ 2 

OT 

2 i!j£_ 


8 

ab 

3 2 W 

3^3 

^ 3X3Y 


"N 

J 


(A. 19) 




j _ 3 2 ~w 

r 3 y 2 

-2^ 


3x3 r 



3 2 w 

3 C 2 

9 2 w 

9r) 2 

3 2 w 

3£3 11 


(A. 20) 


U = ia-Oa-n)^ + i(n? )(i-n)u 2 + |(i+?)(i+n)u 3 

+ |(l-S) (l+C)u^ (A. 21 ) 


V = 1(1-0 (1-1))^ + |(HC )d-n)T 2 + l(nO(i+n)v 3 

+ |(1-0(1+ti)v 4 

w = N-jW-l + n 2 p^ + n 2 p^ + n 4 w 2 + k 5 p 2 + n 6 p 2 + h 7 w 5 
+ I 8 P 4 + Kgpg + K 10 W 4 + H ll p l + H 12 P 2 


(A. 21) 
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ry 


Prom the above , ~ , J-J can be found out and can 

3 5 311 95 3 *) 

be written as below: 


/ ^ 

Hi 

j 

ax 


au 


ax 
\ w • / 



Similarly 


k^ k (1+,l) 

k (1 '«) - k (1+5) k (1+e) 


2a 


d+n) 


i-(i-r) 
2b > l : 

H 



k | 

J 

U 2 

i 

V 

I 

u 4 

1,1 ,* 


(A. 22) 


; - i 

an 


a x 



r = 

a v 


,3YJ 



Ijd-n) fed-n) |;d + n) 


2b(!“C) ' 2b^ 1+ ^ 2b^ 1+ ^ 


- k (1+Tl) 


2b 


(i-c) 



( ^ 


T i) 

.■ ■■ 

. ; ? 

' V 2 1 


<j 

1 

V 4 1 
/ 


(A- 23) 


Prom the above 


{£} 


f an 
ax 

a v 
ax 

a u , 3 v 

ST~ + a"x 


2 w 

ax a? 


(A. 24) 


[B ] gL { 6 } 
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' M' 

3 X 

1 

- ■ ' 

f 

2 

a 

■ 

*> 

3 w 

35 

j 

3 w 

2 

3 w 

3 y 

> 


b 

3 T) 


3 2 w 

35 2 


a 2 U, 


_2 3 _5 

a 9- 5 

_2 Li 

■bat] 


2 HI 


8, 3' % 

2 S _^2 

b 3 T) 


2 

a 


3 K 


12 


35 

ai. 


2 

b arj 


3N 

2 __1£ 
a 35 

2 !!l2 

b 3 T) 


/«/' 


(A. 25) 


3 2 N 


35 


2 1 


+ 


9£ 

3 2 N 




92 ", 1 9 " 4 

^Po + — ^ w 2 + 


# N. 


3 5 5 


6 n 2 , 

Po + 


3 2 N 


+ 


3 5 


10 

— w„ 
2 4 


3 5 
3 2 N 
35 
3 2 N. 


2 K 2 


35 
a 2 ! 


-2 6 
2 P 1 


W, + 
2 2 


35 


8 p 3 1 

r p i + 


35 

3 2 N 


35 


2 ^ 


35 5 


- 


3 2 n 

H 4 

3 5 2 2 


3 2 w 


3r) : 


3 2 N. 


1 


3n 


W-, 4- 


2 "1 


>rii ?s‘ 1 

2 p l + 2 pi 


3 n 2 1 


•n 

3 2 N 


3 r] 


2 K 2 


3n 2 

3 2 N„ 


3 2 W , 3 2 1H o 

4 5 o a , d p2 

1 + 1 - Pt + 7 P2 


3 T) 
3 2 ^ 


2 


+ 


3 ~ 2 


Lw, 


ri 


3 ri 


- P 3 + 

2 1 


3n 

3 2 N 


a fl3 


3 r. 


2 2 
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fl 2 W_ 

6>«,an 


3 2 N, 


353n 


+ 


3 2 EL 3% 

2 1 i 5 

^ + h 


3£3ri 
3 2 N 
3?3n 
3 2 U, 


w 2 + 


1 


3£3n 

3 2 BL 

2 

3?3n 
3 2 3ST 


3 2 N 


fi£ + 


3?3ti 


6 ft 2 
" P 2 


3?3n 

3 2 N. 


w 3 + 


§ P 3 + 

3?3n 1 


3 2 IT 


3£3n 




10 


3 2 N. 


3C3n 


w. 


11 ft 4 , 
p i + 


3 2 M 


3?3n 


12 r 4 

Po 


3C3n 


Prom the above 


{*} = 


can be found out, or. 


r _ 3 2 w 


r _ 4_ 3 2 w^ 

3X 2 


” a 2 3C 2 

_ iiss r 

= i 

- 4 3 2 w 

1 3Y 2 ! 


b 2 3 n ! 

„3 2 w 


8 3 2 w 

i -2 

1 9X9Y 

i 

ab 3^a n 

^ J 1 

L J 


{K> = fe] B {6> e 


(A. 27 ) 


The inplane and bending stiffness matrix can be ca-lculated 
by using numerical integration technique. 

The original double integral in cartesian coordinates 


has been transformed to 
1 

/ f ( £,T))d E&T). 
-1 -1 


f 
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n n 


The above integral is equal to Y y H.H.f (a. ,b. ), 

i=l j 1 0 i 

where H. and H. are the weighing coefficients, f(a., b.) 
0 1 i 

are 'the values of the function at specified points 

(a,,b.) and n the number of Gauss points used. The 
3 l 

present integration has been carried out for 6 Gaussian 
points . 
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APPENDIX B 

COMPUTER INTEGRATION OF CONGRUENTLY TRANSFORMED MATRICES ^ 

In the derivation of finite element expressions for 
stiffness or for stability matrices, certain integrals must 
be evaluated. Usually one of two following methods is 
Employed for this evaluation: 

(a) Analytical evaluation, which involves a considerable 
amount of tedious algebra. 

(b) Numerical calculation, which may introduce truncation 
errors at an early stage of computation. 

A method for the evaluation by the computer of 
analytical expressions for the integrals is presented. 

It avoids the tedious algebraic work, and it evaluates the 
exact analytical expressions for the integral. A FORTRAN IV 
program is written for the evaluation of these integrals for 
the stiffness of the element. 

In structural analysis by finite elements usually 

either the displacement field or the equilibrium field is 

represented by a polynomial. This usually gives rise in the 

derivations to the following type of integral which generates 

( 0 V 

the element of a matrix A K 

A^ = / C B d R • (B.l) 

R 


+ A.P.Kabaila and P. Beckers (81). 
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The terms of matrix B consist of polynomial terms, 
matric C is of constant terms and R indicates the domain 
of integration. Sometimes under the integral one also has 
factors X,Y or Z or all. 


In plate bending, plane stress or stability inte- 
grals, eq uation (C.l) becomes 


CO) r r T (Y-v-"*-!*) (Yy“^-) 

A ( ' = J J B C B X S Y 1 dXdY 


(B.2) 


where and y v are equal to or greater than 1, depending 

on whether there are factors of X and Y or not (in addition 

T 

to those contained B C B). 


The method presented is for the evaluation of 
equation (B.2). The method yields theanalytieally 
exact value, except for the round-off errors, and with 
slight modification can be applied to problems involving 
triple and single integration. As an introduction to the 
solution of equation (B.2), consider the follow ing matrix 
operation with matrices of constant coefficients 


A = B 


T 


C B 


(B.3) 


A typical element of matrix A can be written as' follows 



I I 


k SL 


ki °k£ 



(B.4) 


where letters with subscripts are the elements of 
corresponding matrices. Equation (B.4) follows simply from 
the rules of matrix algebra, because 
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and 


(0 b > kj = f °W b * j 


a ij - l (b \k (Cb V 


l • V «'i | C kJ. hj 


^ | b ki ^k£ \ j 


q. e. d , 


Expression (B.4) provides an algorithm which is very 
simple to programme, though it is not necessarily an efficient 
algorithm as far as computation time is concerned. 


Equation (B.4) can be modified for the evaluation 
of equation (B.2), since a typical term of matrix B of 
equation (B.2), denoted by B. can be written as follows: 


d . . g. . 

B-. = b.. x ^ 

ID ID 


(B.5) 


where b is a constant factor, generally a floating 

-L. J 

point number, and d. . and g.. are exponents of the two 

coordinates, usually fixed point numbers. The values of 

b- d. . and g. . are the data for the purpos es of this 
ID ID ID 

discussion. Substituting B. . in place of b. . in 

J tl 

T 

equation (B.4) one can write the elements of B C B of 


equation (B.2) as follows: 


T 

(B C B) ± j - l l ~ ki 


l l CL 


X 


(^,+doJ (gv,+g„J 


ki 


^ Y ki *3 (B. 6 ) 



no 


Substituting equation (B.6) into (B.2) and observing 
that the summation can be done before integration, one 

obtains 


_ (0) _ y y c v b 

a ij " J. | C k£ \i R 


(d v . +d„. +v -1) (g. . +g. .+ Y -l) 

. [j t kl x 7 kl y 


dxdy 
(B.7) 


where a {j ^ are ike elements of matrix of equation 

(B.2). Equation (B.7) is the basis of the algorithm. 


The integration 


I = // x^ n_1 ^ y (m-i ) dx dy 
R 


(B.8) 


is easily carried out analytically where 


n = d ki + + r x 


m = Ski + g + ^ 


(B.9) 



Figure B.l: Example of an integration domain. 

For instance, if the domain of integration is a triangle 
in skew coordinates (Figure B.l), then 
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J J y dX dY = F(n,m) h n S m 


(B.10) 


where 


F(n,m) = f (“) Jhil— = y (“) -tliilL (B . 11 ) 

3-0 3 m(n+j) i=0 1 n(m+i) 


and 


m(m-l) . . . . (m - j + l) 

(?) 

J .i! 


Substitution of equation (B.10) into (B.7) yields) 

( 0 ) 


= l In. b v.- b .^ F(n,m) h n 2^ 
±3 k 


= “ ki M 


(B . 12 ) 


Concluding, it is seen from the example that the 

above greatly facilitates computer programming for the 

stiffness and some other matrices of finite elements. It 

also allows a i substantial saving to be made in the 

analytical work associated with the programme. Furthermore, 

the. likelihood of mistakes is greatly reduced. The co-mputer 

time can be greater when using this method as compared 

with a subroutine with algebraically evaluated and, term-by- 

term programmed, elements of . This, however, is a 

small penalty (if any) when one is developing and testing 

(0) 

new elements. It should be stressed that the values of 
evaluated by the method are exact, except for the rounding 
off errors. The truncation error, which is often encountered 
in numerical integration, is nonexistent in the present method. 
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S TIFFNESS MATRIX FOR STRAIGHT B EAM WITH RESPECT TO 

ECCENTRIC AXES 

In this appendix, the stiffness matrices for a 
straight beam parallel to X-axis and parallel to Y-axis 
respectively with respect to eccentric axes are given in 
tabular form. In addition to the assumption of a homogeneous 
isotropic elastic material, assumptions are made which are 
of the same order of accuracy as those made for the shell. 

They are: 

(1) Shear deformation is neglected (the 1 Kirch*. off- 
Lcve Approximation’). 

(2) The beam cross-section possesses two axes of 
sy mm etry so that the shear centre coincides 
with the intersection of the principal 
axes. 

Figures C.l and C.2 show the intersection of the 
respective edge beams with the shell. The directions of 
the global axis system are clearly shown in each figure. 

Tables C.l and C.2 give the stiffness matrices of 
straight edge beam parallel to X-axis and Y-axis respec 
tively with respect to eccentric axes. 
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APPENDIX D 


S KEWED BOUNDARY COEDITIONS 

In order that the skewed conditions may be treated 
as normal constraints, the coordinates of those nodes 
where skewed conditions are specified, are transformed. 
This transformation procedure is analogous to the trans- 
formation from local to global coordinates. 


X* 



Figure D.l: Skewed Kinematic Constraint. 

A prime indicates the skewed coordinate system. 
The transformation for displacements at the ith node as 
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[S jL ] is a simple point transformation involving the 
direction cosines which relate the global and skewed 
systems. The transformation for the entire nodal displace- 
ment vector as 

{r } = [S] { r> } (D.2) 

Here the matrix [s] is of the form 



where [i] is an identity matrix of the same order as [sj 
and the number of submatrices on the diagonal of [s] is 
equal to the number of nodes in the ass emblage. It 
follows that the order of CS ± ) and [i] is equal to the 
number of displacements degrees of freedom at each node. 
Note that there is one [S ± ] matrix in [S] for each node 
with skewed constraints. The resulting transformations for 
the stiffness and loads are 
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f>] = [S] T [K] [S] 
{R'} = [S] T {R } 


(D.3) 


Note that [K]' like [K] is symmetric. 

The procedure to transform skewed boundaries is 
often performed before the individual element stiffness 

and loads are assembled. 



